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Figure 1: General sampled data system 


1 MARSYAS Overview 

MARSYAS is a computer-aided control system design package for the simulation and analysis of open loop 
and closed loop dynamic systems. Prior to the summer of 1991, MARSYAS functios provided simulation 
of sampled data (mixed discrete-time and continuous-time) systems; however, sampled-data system analysis 
of stability, frequency response and transfer function coefficients was not implemented. Analysis of purely 
continuous-time systems was available, but the underlying computational procedures did not employ several 
recent advances in numerical techniques. 

This report outlines the numerical and theoretical basis behind the MARSYAS functions developed 
during the summer of 1991. In particular, the numerical computation of the matrix exponential e A = 
+ 88 described i n [11] and the numerical computation of the finite system zeros of 

a dynamic system (continuous- time or discrete-time) as discussed in [4], [6], and [8] are presented. These 
two numerical functions and their associated numerical tests comprise the bulk of the work done by the 
author while working in tandem with John Tiller (BCSS) and D. Pat Vallely (NASA-MSFC). The analysis 
of sampled-data systems is discussed in Section 2; the numerical computation of system zeros is discussed in 
Section 3. 


2 Discretization of LTI Continuous Systems 

A general sampled-data system is shown in Figure 1. The system has an external continuous-time input 
r c (£) and an external discrete-time input r^(/), and maps these signals to continuous-time outputs y c (t) 
and to discrete-time ouputs u c (£). ^ ue to the use °f A/D and D/A converters, sampled-data systems are 
nonlinear and time- varying in general. That is, there is no transfer function from r c (t) to y c (t), even if the 
constituent subsystems are linear and time-invariant (LTI)! However if the subsystems are LTI, the behavior 
of the A/D and D/A converters (u c (f) = Ud(kT), t € [kT, JcT + 7 1 )) implies that one may obtain discrete- 
time transfer functions Yd(z)/ Rd(z) and Ud(z)/Rd(z), This is accomplished by computing matrices F and 
G such that a continuous time system x ~ Ax -f Bu may be equivalently represented at sampling times 
0, T, . . . , k T, ... as x(k T + T) = Fx(kT) 4* Gu(kT). This conversion is obtained via the relations F = e AT 

, G — ^ j B ) where the matrix exponential e AT is defined as e AT = 7 + AT -f (AT) 7 / 2! -f * * % 

Van Loan [10] observes that F and G may be simultaneously computed as 

that is, only a single matrix exponential is required. 

Computation of the matrix exponential is not a trivial task. (See, for example, [7].) However, Ward 
[11] has proposed the use of Pade approximations with preconditioning for this purpose. The algorithm is 
shown in Figure 2. This procedure was tested at MSFC during the summer of 1991 in several numerical 
experiments; computed eigenvalues of F = e A were compared with the exponentials of the eigenvalues of 
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Figure 2: Computation of the matrix exponential 

A y with good agreement. Similarly, matrices A = V were constructed so that computed values of F 
could be compared with exact (known) solution values, again with good agreement. Experiment models 
had up to 109 states with wide variation in coefficients. Little difference was found between using 3 rd order 
Pade approximations and eigth order Pade approximations; hence, for the sake of computational speed, 
MARSYAS presently uses 3 rd order approximations. 


3 Facto red Transfer Fun ction s 

Numerical experiments involving discretized sampled-data systems with fast modes revealed that existing 
MARSYAS software was not adequate for the computation of associated system zeros. Prior to the summer 
of 1991, MARSYAS computed the zeros of H(s) = C(sl — A) 1 B + D as 
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Unfortunately, discretized systems obtained from continuous time systems with fast modes (relative to the 

f —A —B 1 . 

sampling interval) are nearly uncontrollable; that is, the matrix CD i S near ^ s * n 6 u ^ ar - 

EISPACK routine balanc [9] and Ward’s GEP balancing procedure [12] attempt to reduce the dynamic 
range of the coefficients in the algebraic eigenvalue problem (AEP) and the generalized eigenvalue problem 
(GEP), respectively, without introducing any roundoff error. Both of these procedures may be divided into 
two steps: (1) permutation and (2) scaling. 

For example, the permutation step gep.perm of Ward’s balancing procedure computes permutation ma- 
trices Pi and P2 such that the transformed matrix pencil 

(M - \N) = Pi(M - *N)P 2 
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may be conformably partitioned as M — 0 A/22 A/23 and N — 0 N22 N22 where A/n, 

[ 0 0 m 33 J L 0 0 % J 

Yu, A/33 and /V33 are upper triangular and A/22 and N22 form a reduced order permutation-irreducible 
GEP. Following gep_perm, the scaling step gep_scale of Ward’s balancing procedure computes diagonal 
matrices D, = diag(r^’ > , . . . , r*-’), t = 1,2, where r is the machine radix. D x and D 2 are selected so that 
the elements of the matrix products D X MD 2 and D\ND 2 are of approximately the same magnitude. (The 
use of powers of the machine radix r allows these matrix products to be computed without roundoff.) 

The AEP balancing procedure balanc may be similarly divided into a permutation step aep_perm 
and scaling step aep.scale. However, balanc requires that P x = P 2 -1 = P2 and D 2 - D x = 

diag(r _i i > ,...,r _i " ) ). The use of similarity transforms Pi and D\ preserves the N — / structure of 
the AEP. 


Q-2 













Algorithm 

Preconditioning 

Computed Zeros 

ENVD 

none 

-1904, (-1500 ± j'8.744 • lO 4 ), -336.2,0,0 

EISPACK 

none 

-5.468 • 10 s , (2.688 ± >4.720) • 10 & , 
(9.913 ± j'9.359) • 10 7 , -1500, -1500 

ENVD 

zgep.scale 

-1500,-1500 

EISPACK 

gep^perm, gep_scale 

1500 ± j'2.7847 • 10 -5 


Figure 3: Numerical results for Example 3.1. 


Since Ward has developed an effective preconditioner for the generalized eigenvalue problem [12], it 
was decided to modify MARSYAS to use this balancing procedure in tandem with the QZ iteration [8] as 
implemented in EISPACK. This procedure resulted in great improvement in the accuracy of MARSYAS 
computed results. Ward’s balancing procedure may be applied without roundoff error and will usually 
provide improvement in the computed results. However, the QZ iteration fails to isolate zeros at infinity, 
and so these zeros may be perturbed in the Riemann sphere to large (but finite) zeros of arbitrary phase. 

Emami-Naeini and Van Dooren [4] propose a procedure for the solution of the zero-computation gener- 
alized eigenvalue problem that reduces the original problem to one of lower degree that has the same finite 
zeros as the original system but with no zeros at infinity. This code has been tested on several industrial 
system models with great success. However, implementation of this algorithm (as available from netlib 
at Oak Ridge National Labs) proved disastrous when applied to a model of the oxidzer-preburner valve 
dynamics of the SSME. This failure was due to widely varying magnitudes in the coefficients in the system 
model. As a part of the 1991 summer faculty program, a new balancing procedure for the zero-computation 
generalized eigenvalue problem 
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was developed as a preconditioner to the Emami-Naeini /Van Dooren algorithm; see [6] for details. Use of 
this procedure in tandem with modified procedures from [12] and EISPACK provides improved numerical 
robustness as shown in the following example. 


Example 3.1 Consider the following seventh order single-input, single-output dynamic system based upon 
the space shuttle main engine oxidizer-preburner valve dynamics. The system coefficients are an = —7000, 
ai 2 = —2.5 x 10 7 , an = —1.82943 x 10 9 , a 2 i = 1, 032 = 6.4 x 10 5 , 033 = —2,240, 034 = — 032 , a 43 — 1, 
a 54 = 9.03934, a 65 = 225,449, a 66 = -3,000,a 67 = -2.25 x 10 6 , a 76 = 1, b u = 1.44813 x 10 8 , c 15 = 
1.26582, and all other matrix entries are zero. The EISPACK implementation of the QZ algorithm and the 
Emami-Naeini/ Van Dooren (ENVD) algorithm were applied to this system with either no preconditioning, 
Ward’s balancing procedure gep_perm, gep_scale, or the zero- computation generalized eigenvalue balancing 
procedure zgep.scale. A summary of the numerical results is in Figure 3. (Ward’s balancing procedure 
[12] correctly isolates six generalized eigenvalues at infinity; the correct finite system zeros are s = 1500 ± 
j2.7847 x 10 -5 . The zgep_scale scaled system coefficients are an = —7000, ai 2 = —762.9, a 17 = —1744.7, 
a 21 = 32,768, a 32 = 78.125, a 33 = -2240, a 34 = -156.25, a 43 = 4096, a 54 = 72.31472, a 65 = 27.52, 
a 66 = -3000, a 67 = -34.33, a 76 = 65,536, b u = 1.0789 and c i5 = 0.6329. 
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